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We show theoretically that in the external magnetic field like charges on top of graphene monolayer 
may be mutually attracted to form macro-molecules. For this to happen graphene needs to be 
in Quantum Hall plateau state with local chemical potential being between the Landau levels. 
Graphene electron (s) gets localized in the middle between charges and provides overscreening of 
Coulomb repulsion between the charges. The size of the resulting macro-molecules is of the order 
of the magnetic length (~ 10 nm for magnetic field 10 T). The possible stable macro-molecules that 
unit charges can form on graphene in magnetic field are classified. The binding survives significant 
temperatures, exceeding mobility barriers for many ionically bond impurities. The influence of 
possible lattice-scale effects of valley-mixing are discussed. 

Tuning the doping of graphene or the magnetic field, the binding of impurities can be turned 
on and off and the macro-molecule size may be tuned. This opens the perspective to nanoscopic 
manipulation of ions on graphene by using magnetic field and gating. 
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I. INTRODUCTION 

More than ten years since discovery of graphen^, 
the first genuinely two-dimensional material, are 
marked with huge body of research. Its applica¬ 
tions range from cancer-drug deliverjl^lSI to novel 
electronic devices. Still, many of proposed uses of 
graphene depend crucially on its interface interac¬ 
tions with other compounds and impurities. Typi¬ 
cally, molecular dynamics and DFT methodP"^ are 
used to model the interaction of graphene with sub¬ 
strate and impurities. It was shown that there are 
two types of impurity bonding to graphen^ cova¬ 
lently bond impurities (e.g. H, CH3, F, OH, O) and 
ionically bond impurities (e.g. Na, K, Cs, Cl, Br, I). 
Covalently bond impurities act similar to defects in 
graphene structure and typically have large mobil¬ 
ity barriers, while ionically bond impurities have low 
mobility barriers (typically lower than room temper¬ 
ature) and act as mobile electric charges put on top 
of graphen^^ESmi, 

Interaction between impurities has attracted sig¬ 
nificant experimental and theoretical interest. It was 
showrfi^ that short-range covalently bond impuri¬ 
ties tend to form bound clusters due to a fermionic 
Casimir effect with binding energy comparable to 
room temperature scale at distances below 2 nm. 
Related results were presented in ref. 0 where di¬ 
lute adatoms were shown to have tendency towards a 
spatially correlated state with a hidden Kekule mo¬ 
saic order. On the other hand, this effect cannot 
compete with Coulomb repulsion of like charges if 
the impurities are charged, as is the case for ioni¬ 
cally bond impurities. 

Graphene is known to be exceptionally suscepti¬ 
ble to magnetic field with divergent diamagnetic sus¬ 
ceptibility at the Dirac point and strong non-linear 


effects in magnetizatioifi^RSl. Magnetic field creates 
large energy gaps near the Dirac point that reduce 
intrinsic screening and help to localize electrons in 
the external potentials. 

Without magnetic field, graphene will p roduce 
non-linear screening of external chargea^^R^. When 
magnetic field is applie d, th e ordinary screening is 
significantly suppresse d^^bQ I and electrons become 
localized. 

In this work we discuss the interaction between 
charged impurities near the surface of graphene in 
the strong magnetic field and show that in some 
cases the binding due to electrons in graphene can 
compete with Coulomb repulsion and stable nano¬ 
molecules may form. The motivation for this work 
comes not only from ionically bond impurities men¬ 
tioned above: Epitaxial graphene on Si-terminated 
SiCP^HMl gg important feature of positively- 
charged donors appearing dynamically in the “dead” 
carbon layer just below graphend^^H^Sl This happens 
due the dominant effect of quantum capacitanc^^ 
and due to strong changes in the density of states 
near the Dira c point under the influence of mag¬ 
netic fielcP^^^. The process looks as the appearance 
of localized holes in the insulating “dead” layer due 
to electrons transferred to graphene. Consequently, 
the charge transfer might dynamically create extra 
localized states and extra electrons to keep the sys¬ 
tem in a robust v 2 quantum Hall plateau state. 
Since the charge transfer is reversible, the holes be¬ 
low graphene can be considered on the same foot¬ 
ing as dynamical charged impurities. The process 
of charge transfer allows us to talk about the well- 
defined local value of the electron chemical poten¬ 
tial. The chemical potential is chosen to be in the 
gap between the undisturbed Landau levels. 

Another possible context is to consider the 
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well- know n semi-classical picture of Quantum Hall 
effecPI^ (QHE). It assumes smooth (on the scale 
of magnetic length) external potentials and derives 
the presence of compressible regions where electrons 
screen everything and the incompressible regions 
that have exactly integer local filling factors. For 
non-smooth potentials created by impurities close 
to graphene sheet, the semi-classical approach is not 
the full truth. The self-consistent semi-classical ap¬ 
proach gives a smoothly-varying effective potential 
that could be locally treated as a chemical potential 
for our purpose. Our results imply that there is still 
some life in the incompressible region due to screen¬ 
ing and possible over-screening of point charges, and 
thus the local filling factor of incompressible regions 
could slightly deviate from integer values. 


For definiteness, we consider the positively 
charged ions below, but exactly the same description 
applies to the negative ions if one makes a particle- 
hole transformation and reverts the sign of chemical 
potential. 


We show that when several charged impurities 
are at th e distan ce of the order of magnetic length 
Ib = y/fi/ (eB), they can form stable molecules 
bound by the electrons in graphene. Depending on 
its local chemical potential (taken in the region of 
Quantum Hall plateu), graphene may bind either 
positive ions, or negative, or produce no significant 
binding at all, the results are summarized in Table 
and Figj^ The described effect may be called either 
as over-screening of Coulomb repulsion of impurities 
or as long range covalent bonding of charged impuri¬ 
ties. A related “overscreening” effect exists for the 
charging of a quantum dot!^ in colloidal systems a 
similar effect is called ’’charge inversion”!^. 


Graphene in magnetic field hosts the electron hy¬ 
bridization cloud that can lead to attraction exceed¬ 
ing the Coulomb repulsion of same-charge ions, see 
Figg Significant binding of charges occurs only 
when the electron cloud is centered in between the 
charges (r 3 on Figg. Simil arly to the signatures 
of Quantum Hall Effect (QHE)P^ we show that the 
effect may survive the room temperature for mag¬ 
netic fields of order 10 T. 


The paper has the following structure. In section 
[Tl| we introduce the formalism and present the re¬ 
sults for bound states bound with 1 electron; in sec¬ 
tion m the results are generalized to multi-electron 
bound states; in section]^ the obtained results are 
combined in a phase diagram; in sectionthe qual¬ 
itative effects of potentials that are non-smooth on 
the lattice scale are discussed; conclusions are given 
in section ED 


II. GRAPHENE WITH CHARGED 
IMPURITIES IN MAGNETIC FIELD 


It is well-known that the one-particle energy levels 
of an ideal graphene in ma gnet ic field B are given 
by degenerate Landau levelJlElHM 


En = sign{n)EB^/2\n\ , n gZ 

( 1 ) 

where vp ~ 10®m/s, 



Ib = \/hl{eB) 

« 26/a/H/(T esla) nm 

(2) 

Ib 

:! 26v/H/(Tesla) meV. 

( 3 ) 


Each Landau level (LL) is degenerate with density 
area. Working with finite area, it is 
convenient to use the basis of Landau wave-functions 
in polar coordinates. Defining an oscillator radial 
eigenfunction: 


1 -Fi [-n] \m\ + 1; 4) 

we have for the wave-functions { m < n , n > 0 ): 




( 4 ) 




ig ™+|m| (r) 


The solutions near the second AT-point are ob- 
tained by acting with ia 2 on the spinor: I ^ I — 


When the Coulomb impurity potential is present, 
the orbital (index to) degenerac y of Lan dau levels is 
liftecP^. This has been calculate d^^F^ps i demon¬ 
strated experimentalljff^ for one Coulomb impurity. 
Below we consider several impurities. 

Consider a superposition of Coulomb potentials of 
the form 


Uif) = - 


E 


47reoe -y/(r — fi)"^ -\- d? 


( 5 ) 


Here the parameter d is a vertical displacement of 
impurity from the graphene sheet, but it can also be 
used to model a finite localization length of impurity 
wave-functioiP^. Here, there are two physical cases 
that may be considered: 

• The wave function of the localized impurity 
exceeds the graphene lattice scale, d 3> a, but 
still d ^ Ib- This will be the main topic 
of the present work. This may be relevant 
to the case when localization of impurity is 
determined by Moire pattern formed by in¬ 
commensurate substrate and graphene. For 
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FIG. 1. Electron densities for the lowest sublevels of 0-th and 1-st Landau levels in the field of 3 ions forming 
equilateral triangle with sides r. Essential binding occurs for r < 31 b where the electron wave-function forms a single 
lump in the center. 
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EIG. 2. Plot of 1-particle energy sub-levels for sev¬ 
eral Landau levels as a function of distance between two 
positively charged Coulomb centers. Here q = 0.4 and 
d — 0.05. A remarkable (but obvious) fact is that if all 
of the sublevels are populated, the distance-dependence 
of the total energy disappears. 



EIG. 3. Zero-temperature phase energy plot for dilute 
mobile charged impurities. Phase with lowest AEwn is 
favored. Legend corresponds to Table |I] The chemical 
potential should be understood as taken relative to Lan¬ 
dau level edges, which in this plot are renormalized to 
free-particle values of 0 and %/2 here. 


epitaxial graphene on Si-terminated SiC this 
scale is estimated as d ~ 2 nm. Then the 
impurity potential is smooth on the lattice 
scale (note that it is still sharp on the scale 
of magnetic length) and single-valley contin¬ 
uum approximation works well, we may ne¬ 


glect the inter-valley scattering. This means 
that the valley-mixing splitting is less than 
interaction-induced energy and multi-electron 
wave functions are determined by long-range 
Coulomb interactions. This case is most uni¬ 
versal since it depends only on the magnetic 
length scale and effective coupling constant. 
Additional universality comes from the fact 
that the dependence of binding force on the 
effective range d of impurity wave function is 
weak, see Fig|^ 

• If the impurity is more like a point charge, 
then the impurity potential is sharp and 
non-universal lattice-scale details come into 
play and cause the valley mixing. Continu¬ 
ous approximation Eq.Q breaks down near 
the impurity center since sharp potentials 
cause inter-valley scattering. This is dis¬ 
cussed in section 0 We discuss several lattice- 
related effects, but the main conclusions of 
the continuum-mo del considerations still hold 
true. 


The equations for single-electron energy levels in 
graphene in the magnetic field B and any Coulomb 
potentials can be rewritterP^ in units of magnetic 
length Ib, magnetic energy Eb and dimensionless 
coupling 

a = / {A'KtQhvFC-es) = 2.19/eeff (6) 

with the effective dielectric constant CeS = (ei-|-e 2)/2 
coming from substrates on both sides of graphene 
and from graphene by itself. For example, a « 1 on 
Si02 substrat^i^ and a « 0.4 on SiC. In dimension¬ 
less units the equation to solve 


[(Ji{idx-y)-i(T2dy\^ 


E-j^ , ! 

( 7 ) 
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For “molecules” bound by more than one electron 
we will account for e-e interactions below. For per¬ 
forming the computations we evaluate the matrix 
elements of the impurity potential in the basis Q 
truncated to several Landau levels (of the order of 
10) and orbital states (of order of 30) and then do an 
exact diagonalization. In the case of smooth poten¬ 
tials it is sufficient to consider only one valley. The 
truncation of orbital states we use corresponds to 
a circular box truncation of space and to avoid un¬ 
physical boundary contributions we had to smoothly 
cut-off the Coulomb potential at large distances (of 
order of 4/^). As a result, for inter-ion distances 
up to 21b, the precision is better than 1% while for 
larger distances the error may get higher. 

The Landau levels that are completely filled do 
not contribute significantly to the energy of ions as 
a function of their separation (but they d o con tribute 
to renormalization of chemical potentiaP^EIl) This 
has been verified numerically and is seen analytically 
in the leading order of perturbation theory in the 
potential: 

A(ri 2 ) ~ ri) + V{r- r2)|V'i°L) 

m 

= 2^1 = const ( 8 ) 

m 

where we used that the full degenerate set of wave- 
functions, corresponding to a given Landau level 
(enumerated by m) maps to itself under translations 
(up to a unitary transformation), and so the above 
sum is independent of the impurity positions. 

The situation changes drastically if only one or 
several lowest Landau sub-levels are filled with elec¬ 
trons. This happens when graphene is in Quantum 
Hall plateau state corresponding to the chemical po¬ 
tential being in the gap below the band edge (i.e. 
between the unperturbed LLs in 1-particle picture). 

Let us start with two positive charges N = 2. The 
single-particle energy levels of electrons in graphene 
are presented on FiglU Note that only the lowest 
sub-levels are fully meaningful as Filing more levels 
requires to account for the e-e interactions. 

When the distance between the two ions is of the 
order of magnetic length Ib, the lowest energy elec¬ 
tron wave-function is centered in the middle between 
the ions and plays the role of hybridization cloud 
that binds them. Fig. When only this lowest- 
energy state is filled, the strong dependence of en¬ 
ergy on the distance ri 2 between ions appears. Fig. 

creating the attractive force, Fig. If the dis¬ 
tance between charges exceeds roughly 31 b , the low¬ 
est energy electron wave-function becomes centered 
near each of the individual charges and its energy 
depends on ri 2 as E ~ — I/ri 2 . 

In the leading order of perturbation theory, the 
binding force is proportional to a. Notably, the mu¬ 
tual Coulomb repulsion of ions is also proportional 


to a: 

Oi 

T/Coulomb — (9) 

ri2 

Thus, in the leading approximation the distance 
where the attraction would balance the repulsion 
is independent of a. We stress here that when ex¬ 
pressed in magnetic length and energy units, there 
are essentially no free parameters in the problem and 
the Coulomb repulsion of ions is of the same order 
of magnitude as hybridization attraction.!^ It is a- 
priori not at all clear if stable bound states of ions 
can form. Moreover, strong bound states form only 
when the chemical potential /r is near the 0-th and 
I-st Landau levels. When /x is near the other Landau 
levels (see Fig. [^, the ri 2 dependence of electron en¬ 
ergy is substantially weaker leading to weaker bind¬ 
ing of ions. The main problem with bound states in 
the higher plateau states is that they will not survive 
graphene rippling and temperature, thus we concen¬ 
trate on the most robust v « ±2 plateau states be¬ 
low. 

Let us first study the dependence on the coupling 
constant a, see Fig. io. We observe that the 0-th LL 
is almost protected from non-linearity, while a signif¬ 
icant non-linearity appears in the higher LL already 
at a ~ O.I. The 1-st and -1-st levels must have the 
same energies in the leading order of perturbation 
theory as their wave-functions differ only by a rela¬ 
tive sign of sublattice components, but we see that 
these levels split already from a ~ 0.1. This could 
be expected due to large Coulomb field near an in¬ 
dividual impurity. Similar asymmetry is present in 
the results of reffTHl The dependence on the dis¬ 
tance d of the ion from graphene (or, on the local¬ 
ization length of charge wave-function) is weak when 
this distance is much less than magnetic length, see 
FiglJ. 

Below we present a particular example of a = 0.4, 
d = 0.051b , the results are universal for 0-th LL, 
while for the 1-st LL the binding is seen to grow a 
bit faster than linear in a. 

The results for two positive ions bound by one 
electron are presented on Figj^. We observe an ab¬ 
solute minimum in the energy at a distance ri 2 min ~ 
1.51b near the 1-st LL or 1.71b near the 0-th LL. 

If the chemical potential decreases below a criti¬ 
cal value, determined as the energy per electron in 
the bound state, see Tablethere will be no bound 
states for positive ions, but, at some point, we start 
getting bound states for negatively charged impuri¬ 
ties bound by holes. The calculations and results for 
negatively charged impurities are exactly the same 
and obtained by going to hole picture. 

Analogous calculation shows that the configura¬ 
tions with three symmetrically positioned impuri¬ 
ties bound by one electron is unstable: despite an 
appearance of local energy minimum, the configura¬ 
tion would gain energy if deformed to a bound pair 
with the third ion repelled to infinity. 
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FIG. 4. Binding force of two ions bound by one electron: a) as a function of inter-ion distance ri 2 for a — 0.4 and 
d = 0.05Zs; b) as a function of coupling a for d = 0.051 b and ri 2 = 1.51b', c) as a function of d (distance of ions from 
the graphene sheet) for a = 0.4, ri 2 = 1.5Zs 




FIG. 5. Energy (electron energy -|- Goulomb repulsion of ions) of stable macro-molecules. Electron energy is 
counted from Landau level edge (fixed to free-particle values 0 and \p2 for 0-th and 1-st LLs respectively). Here 
a = 0.4 and d = 0.051b. Dotted lines mark the r —>■ oo asymptotic to show stability, a) Two ions bound by one 
electron in graphene; b) Two ions bound by two electrons in symmetric orbital state; c) Three ions forming equilateral 
triangle bound by two electrons in symmetric state 


III. MULTI-ELECTRON STATES 

Now consider multi-electron bound states. The 
two-electron bound state can be in a symmetric 
or antisymmetric orbital state. In a conventional 
molecule this would correspond to spin singlet and 
triplet respectively, but in gr aphene one has an addi¬ 
tional valley degenerac} *^"^ * "^^ * and the full SU(4) sym¬ 
metry is approximately respected (neglecting the 
Zeeman splitting and valley-mixing). Group theory 
tells us that 4 (g) 4 = 10 0 6 in spin-valley space, so, 
there are 10 possibilities to form antisymmetric or¬ 
bital state and 6 for symmetric one. If valley-mixing 
effects discussed in section are important (exceed 
the temperature scale), then orbital-symmetric state 
will be the usual spin-singlet. 

To find a reasonable approximation to the energy, 
we use the variational Hartree-Fock method with the 
basis formed by Slater determinants of low-energy 1- 
particle eigenstates found above. 

Looking at Figj^ we note that near the 1-st LL 
the two lowest single-particle energy sublevels grow 
with ri 2 thus can potentially bind the ions and can 
participate in antisymmetric orbital wave-function. 
For the 0-th LL only the lowest sublevel binds the 
ions while electrons in the higher levels do not tend 
to hybridize, thus the orbital wave-function should 
be symmetric for a stable ion-binding. These con¬ 
clusions were checked by explicit computations and 
comparison of the energy. 

Consider two electrons in the field of a single ion. 
For a = 0.4 we have energies —0.28 , —0.21 per 
electron in 0-th and 1-st LL respectively for sym¬ 


metric state and —0.24 , —0.22 respectively for anti¬ 
symmetric orbital. Thus we expect symmetric state 
near 0-th LL and antisymmetric state near the 1-st 
LL. Calculation shows that one unit-charge ion can¬ 
not hold more than 2 electrons in the lowest-energy 
states: e-e interactions make this too expensive. 

Now we consider two positive charges and two 
electrons. Calculation shows that symmetric or¬ 
bital state is preferred. The resulting hydrogen-like 
molecule is very stable with optimal inter-atomic dis¬ 
tance 1.21b and 1.31b for 1-st and 0-th LL, Fig. [^. 

The same calculation can be repeated for 3 ions 
bound by two electrons, see Figj^. We compare the 
triangle configuration of ions with a linear chain ge¬ 
ometry and find that equilateral triangle geometry 
has lower energy. The described equilateral triangle 
with 2 electrons is prominent for providing the lowest 
possible energy per electron, thus, it is this configu¬ 
ration that appears first in the phase diagram. Table. 
1^ Four ions cannot be bound by two electrons. 

Considering now 3-electron states in Hartree ap¬ 
proximation we found that three electrons cannot 
form a one-centered wave-function to bind any num¬ 
ber of ions since the e-e interaction gets too high and 
the resulting energy gain can by no means compete 
with the energy of far-separated smaller clusters de¬ 
scribed above. The situation is different from the 
ordinary molecules since the shape of wave-function 
is mostly determined by Landau level number and 
not by Coulomb potential of charges. For the low¬ 
est Landau levels that are most robust, the shape 
of wave-functions is essentially circular (similar to 
s-state). The above considerations still do not ex- 
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elude the possibility of larger multi-electron bound 
states with multi-centered electron wave function, 
but quantitative study of such configurations is be¬ 
yond the scope of this paper. 


TABLE I. Phase diagram for N positively charged ions 
bound by n electrons at a = 0.4, d = 0.051b', for 
negatively-charged ions one has to replace /i —>■ —fj,. 
Umin. gives the minimal chemical potential for existence 
of a given phase. Chemical potentials should be under¬ 
stood as taken relative to Landau level band edges, which 
are here renormalized to free-particle values 0 and %/2 
for easy match with free-particle levels. Vij is the opti¬ 
mal distance between ions, i.e. the size of the molecule, 
measured in units of magnetic length Is- E is a minimal 
possible (binding) energy per molecule, which is achieved 
when fi tends to the band edge (p = 0 for configurations 
1 ... 5 and /r = \/2 for configurations 6 ... 10). For 
lower jj, the energy to consider is E — nAfj,. The last 
column, fj, interval, indicate the range of chemical po¬ 
tentials where the phase gives the minimal energy per 


charge (see Eq.( 10 


yo 

A^min 

N ions 

n electrons 

rij 

E 

fj, interval 

1 

-0.93 

3 A 

2 symm. 

1.5 

-1.06 

[-0.93,-0.75] 

[-0.55,-0.38] 

2 

-0.75 

2 

1 

1.7 

-0.52 

[-0.75, -0.55] 

3 

-0.63 

2 

2 symm. 

1.3 

-0.96 

[-0.38 ,-0.08] 

4 

-0.45 

1 

1 

- 

-0.45 

0 

5 

-0.28 

1 

2 symm. 

- 

-0.56 

[-0.08,0] 

6 

0.51 

3 A 

2 symm. 

1.4 

-0.95 

[0.51,0.67] 

[0.95,1.09] 

7 

0.67 

2 

1 

1.5 

-0.47 

[0.67,0.95] 

8 

0.82 

2 

2 symm. 

1.2 

-0.85 

[1.09,1.39] 

9 

1.02 

1 

1 

- 

-0.39 

0 

10 

1.19 

1 

2 anti-s. 

- 

-0.45 

[1.39,1.41] 


IV. PHASE DIAGRAM 

Having studied the simple macro-molecules sepa¬ 
rately, the results may be combined in a phase dia¬ 
gram, see Table|l] In each of the above states one can 
find the energy per electron that bind the molecule 
and thus find a minimal electron chemical poten¬ 
tial /Xmin for such molecule to appear. For a given 
electron chemical potential, several molecule config¬ 
urations may be possible. Let us assume that the 
number of ions and electron chemical potential are 
fixed. For a molecule with JV ions bound by n elec¬ 
trons we compute the free energy gain per ion by the 
formula: 


A 7 ^ _ (I^binding + -^Coulomb) - ftn E- fj.n 

=- N -= 

( 10 ) 

The configurations with minimal AFlion will deliver 
the minimal free energy at zero temperature, see Fig. 
1^ In particular, Fig|^ shows that the state of one 


electron bound by one ion (lines 4 and 9 on the plot) 
is never the lowest-energy state. At high chemical 
potentials (close to the vacuum LL) the state with 

2 electrons per one ion wins (lines 5 and 10). All 
the other states correspond to bound states of 2 or 

3 ions, which are realized for ^ S [—0.93, —0.08] U 
[0.51,1.39] for the example a = 0.4 considered. 

Changing a in the leading approximation just lin¬ 
early scales the phase diagram around the vacuum 
Landau level (0 and V2). As is seen from Fig. 
the binding of molecules near v = 2 should addi¬ 
tionally increase with increasing a due to noticeable 
non-linearity. 

To come to physical conclusions, consider a real¬ 
istic example of randomly positioned charges. For 
example, consider 4 charges and the z/ « 2 plateau 
state. Looking at the Table we can find the states 
that deliver the best binding energy per charge: 

• if the chemical potential is /i G [0.51,0.67], 3 
charges may be bound with 2 electrons in the 
triangle (state iN' 26 ) and the fourth is repelled 
to infinity since no electrons can be bound to 
it. In reality, it would be hard to create such 
state since one needs to surpass a signiheant 
Coulomb repulsion gap to form the triatomic 
molecule. 

• if /i S [0.67,1.02], two diatomic ions (two ions 
bound with one electron) would be the pre¬ 
ferred state (again, it will be hard to create). 
The state 1 will continue to exist as well. 

• for p, > 1 . 02 , a single charge can bind elec¬ 
tron, forming a neutral combination. Now it 
is much easier to form molecules since there is 
no long-range Coulomb barrier. In particular, 
the triatomic configuration (JV^O in the table) 
is now much easier to form, note that this con¬ 
figuration has unit total positive charge. One 
remaining charge will have one electron bound 
to it and this neutral combination will be at¬ 
tracted by the induced dipole moment to any 
charged object. Thus, a 4-atomic molecule 
may be formed. The precise determination of 
its energy is beyond the scope of this paper. 

• for fi > 1.09 the diatomic charged configura¬ 
tions (N 26 ) would shrink in size (from 1.5 to 
1 . 2 ) to form the neutral diatomic molecules 
(N28). 

Analyzing the above example we may draw the fol¬ 
lowing conclusions. If there are randomly-positioned 
same-charge mobile impurities, the active recombi¬ 
nation process would start only when neutral states 
of one impurity and one electron (JV^d or JV29) may 
be formed. For lower chemical potential the stable 
bound states exist but they are hard to form due 
to significant Coulomb repulsion barriers. The most 
relevant recombination channel is to form neutral 
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hydrogen-like molecule (JV^S or JV^S). In general, the 
clustering process is mediated by the induced dipole 
attraction of neutral and charged macro-molecules. 

If we are dealing not with the real ionic impuri¬ 
ties, but with the holes dynamically appearing in 
the surface layer of substrate below graphene, the 
Coulomb barriers for formation of bound states are 
much lower and all the states in Table H] are relevant. 
Their appearance will be governed by the local chem¬ 
ical potential formed by other non-mobile charges. 

It should be kept in mind that all these calcula¬ 
tions make sense when concentration of mobile im¬ 
purities is low: riimp ^ Correspondingly, the 
filling factor of electrons that bind these impurities 
must be very close to the perfect values ±2. These 
conditions could easily be achieved locally in the in¬ 
compressible regions appearing dynamically in the 
QHE physics. All the “extra” electrons are expelled 
out to compressible regions. 


V. VALLEY MIXING AND LATTICE 
EFFECTS 

In case of sharp potentials we need to consider 
the doubled basis, involving the Landau level states 
in both valleys. Naively, the matrix element be¬ 
tween Landau states could be computed by discrete 
lattice summation in the circle around the impu¬ 
rity and with continuum integration in the other 
regions. On the other hand, using only Coulomb 
interaction and Dirac bands in the immediate vicin¬ 
ity of the impurity would be an oversimplification 
and the full (e.g.DFT) study including high-energy 
band is needed. Thus, we follow a phenomenologi¬ 
cal approach and model the impurity as smoothed 
Coulomb potential (as studied above) supplemented 
with point-like scatterers sitting on A and (or) B 
sublattices. According to refUJ positively charged 
impurities (Li, Na,K, Cs) tend to be located in 
the center of graphene hexagons (“h-point”) and 
the negative impurities (Cl,Br,I) prefer to stay on 
top of carbon atoms (“t-point”). Thus, it is natu¬ 
ral to expect sublattice-symmetric effective poten¬ 
tial for positive ions (VAi = Vsi) and sublattice- 
asymmetric potential for negative impurities. In the 
latter case, similarly to ref. [12 it is interesting to 
study the difference of same-sublattice and different- 
sublattice locations of impurities. A general point 
potential is parameterized by its sublattice compo¬ 
nents Vaai, Vbbz, Vabi at each of the impurity posi¬ 
tions Xi- These components are naturally combined 
in 2 X 2 matrix V{xi). 

In the leading order of perturbation theory, these 
potentials shift and split the single-valley energy lev¬ 
els discussed above. Consider corrections to energies 
of I-electron bound states in the field of N impuri¬ 
ties. The wave-functions with the same energy Ag in 
the second valley are obtained as k' = icc^^K- To 


discuss the lifting of this degeneracy in the leading 
order of perturbation theory we have the Hamilto¬ 
nian H = ^ \ with matrix ele- 

\ T'*2) ^0 + ^22 J 

ments 

N 

Vn=Y.'SH^j)V{xj)'Six,) (II) 

i=i 

V^2=J2'^H^,)Vix,)za2'fix,) ( 12 ) 

i 

N 

^^22 = '^Hxj)cr2V{Xj)a2'^{Xj) ( 13 ) 
i=i 

The energy eigenvalues are: 


Vii T V22 

E = Eo+ ± 


f Vll-V22 


lU 


(14) 

Note that the sublattice components of the wave- 
functions for the zeroth and the first Landau-levels 
are strongly asymmetric: one sublattice component 
is dominant for a chosen valley, see Fig[^ Hence, 
with a good accuracy, valley matrix elements are 
proportional to sublattice matrix elements: Hn ~ 
Y.^VBB^\'^{xi)\'^, V 22 ~ Li T'A^i|4'(a;i)p and V 12 ~ 
Li r 4 Bi|'I'(a:i)p. For negatively-charged ions situ¬ 
ated on the same sublattice, we get the main non¬ 
zero matrix element Vn = 2V. This shifts down the 
energy in one of the valleys by 2|F|. If impurities are 
on the different sublattices we get Vn = V 22 = V, 
which gives an equal twice smaller energy shift \V\ 
for both valleys. Thus, location of impurities on the 
same sublattice is energetically preferred and the 
corresponding electronic wave-function prefers this 
sublattice (assuming that H < 0 ). 

For positively-charged ionic impurities, one ex¬ 
pects Vii = V 22 due to preferred location of impurity 
in the centers of hexagons. Then the sublattice¬ 
mixing matrix element V 12 plays the leading role 
in lifting the valley degeneracy. The preferable 
wave-function would then be a symmetric or anti¬ 
symmetric combination of both valleys or sublat¬ 
tices. 

Both effects advocated above lead to an extra gain 
in binding energy that would depend on the dis¬ 
tance between the impurities, thus, producing an 
extra force. The distance dependence is governed 
by |'I'(a;i)^|, which, as already stated, is dominated 
by only one sublattice component, see Fig. 

Figurej^illustrates that the effects related to scat¬ 
tering on sharp potentials do further increase the 
binding force. The 2-electron states discussed be¬ 
fore essentially fill the same orbital state with two 
electrons, hence, the discussed correction just dou¬ 
bles. The matrix elements discussed above scale 
with magnetic field as \ilj‘^{xi)\ ~ 1//^ ~ B. The 
resulting extra contribution to energy scales as B 






''12 its) 


FIG. 6. Plot of sublattice components of wave-functions 
for lowest-energy sublevels of O-th and 1-st Landau Lev¬ 
els in the field of two charges at distance ri 2 . Wave- 
functions are evaluated at the location of charged impu¬ 
rities. 

and the force scales as . This is to be compared 
with scaling of the main term in the energy. 


VI. DISCUSSION AND CONCLUSIONS 

To conclude, we have shown that for a signifi¬ 
cant range of chemical potential values inside the 
gap between Landau levels the charged impurities 
(or donor states in the substrate surface) can form 
stable bound states. The optimal distance between 
charges in the bound s tate is of the order of magnetic 
length IB = The binding energy scales as 

- and the binding force at optimal distance 
scales as F ^ Eb/Ib vpeB. 

The above results were obtained for an ideal 
monolayer graphene sheet at zero temperature. The 
realistic graphene may have other non-mobile impu¬ 
rities, ripples, corrugations and finite temperature. 
Clearly, mobile charge impurities can equally well 
form bound states with non-mobile ones. Simulta¬ 
neously, charged impurities may introduce smooth 
inhomogeneities in the chemical potential leading 
to replacement of chemical potential /r with a lo¬ 
cal chemical potential /i -I- Bimpurities in our consid¬ 
erations. It is also important to note that sufficient 
amount of mobile charged impurities leads to screen¬ 
ing of potential landscape thus making it flatter on 
large scales. This may be one of the keys to under¬ 
standing of exceptionally precise Hall quantization 
in epitaxial graphen^^. 

The temperature and short-range impurities lead 
to level-broadening. As is clear from Fig. the 
binding appears when the lowest Landau sub-level 
is filled, while the next ones are empty. The split¬ 
ting b etween thes e levels is of the order of oBb « 
SOOofy/B/ (Tesla) K (see Fig. [^. The splitting of 
levels is twice smaller near the 1-st LL, but in this 


case the second smallest LL is also attractive and 
population of this level does not spoil the binding. 
So, with a « 0.4 our results must survive the room 
temperatures and small amount of short-range im¬ 
purities for magnetic fields above 10 Tesla and even 
higher temperatures at larger fields. These conclu¬ 
sions are also supported by the experiment of ref. 
m- Note that at room temperature (and even be¬ 
low) many types of ionic impurities are mobil^. 

Another effect of finite temperature is an entropic 
contribution coming from the approximate spin and 
valley degeneracy. This effectively decreases the 
energy by fesTlnn, where n = 4 for one-electron 
bound states, n = 6 for two-electron symmetric 
states and n = 10 for two-electron anti-symmetric 
states. These numbers are to be changed if the de¬ 
generacy lifting due to local lattice effects discussed 
in section |V] or Zeeman splitting exceeds the tem¬ 
perature scale. 

An important aspect in graphene is rippling and 
corrugation^23®6| argued in ref. [47], corruga¬ 
tions in graphene may be described by fluctuations 
in perpendicular magnetic field that lead to consider¬ 
able broadening of non-zero LLs. At the same time, 
the 0-th LL is m'otected and mainly broadens due 
to temperature^. Our main results correspond to 
chemical potentials in the gap above or below the 0- 
th LL, thus the effect of corrugations is expected to 
be moderate. At the same time, corrugations may 
kill any weak binding effects that might occur in the 
higher QHE plateaux. 

To summarize, we have shown that graphene in 
magnetic field can mediate strong attraction of like 
charges put near graphene. The resulting size and 
configuration of macro-molecules depend on mag¬ 
netic length and local chemical potential and thus 
can be easily changed by tuning the magnetic field or 
doping of graphene. The results are expected to sur¬ 
vive significant temperatures. This opens a perspec¬ 
tive to nanoscopic manipulation of ions on graphene 
by using macroscopic tools and provides further in¬ 
sight onto the structure of incompressible regions in 
QHE physics beyond the semiclassical approxima¬ 
tion. Apart from that, the results might shed light 
on the microscopic structure of potential landscape 
in SiC epitaxial Quantum Hall Effect devices. 
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